// --------------------------------------------------------------------
//
//  File:        OAk8D3.rtf [Optimal Ate pairing for twisted Edwards curves with quadratic twists]
//  Date:        Oct. 10, 2018
//  Last update: Oct. 10, 2018
//  Description: Magma implementation of the optimal ate pairing on twisted Edwards curves with quadratic twists.
//               Embedding degree k = 8, CM discriminant D = 3
//
//  (c) 2018, Georgios Fotiadis & Chloe Martindale
//                 
//  Georgios Fotiadis, University of the Aegean, Greece, email: gfotiadis@aegean.gr
//  Chloe Martindale, Technische Universiteit Eindhoven, the Netherlands, email: chloemartindale@gmail.com
// --------------------------------------------------------------------
//  Polynomial Family
//  p(x) = (3*x^16 − 9*x^12 + x^10 − 2*x^9 + 16*x^8 − x^6 + 5*x^5 − 13*x^4 + x^2 - 5*x + 7)/3
//  t(x) = x^8 + x^5 − x^4 − x + 2
//  r(x) = x^8 - x^4 + 1 
// --------------------------------------------------------------------
//  Optimal ate pairing: a(Q,P) = [ f_{x,Q}^{x + p} (P) f_{x,[x]Q}(P) h_{[p^2]Q,[xp]Q}(P) ]^((p^k - 1)/q)
//  Here H = h_{[p^2]Q,[xp]Q}(P)
// --------------------------------------------------------------------

//Weierstrass to Montgomery
MontPoint:=function(P, alpha, s)
  x:=s*(P[1] - alpha);
  y:=s*P[2];
  MP:=Vector(2,[x,y]);
  return MP;
end function;

//Montgomery to Edwards
EdPoint:=function(P)
  X:=(P[1]/P[2]);
  Y:=(P[1] - 1)/(P[1] + 1);
  Z:= 1 ;
  T:=(X*Y)/Z;
  EP:=Vector(4,[X,Y,Z,T]);
  return EP;
end function;

//Computes the components in the value H of the optimal ate pairing
hh:=function(Q1, Q2, P, a)
  x1:=Q1[1]; y1:=Q1[2];
  x2:=Q2[1]; y2:=Q2[2];
  xP:=P[1]; yP:=P[2];
  if x1 eq x2 then 
    h:=xP - x1;
  end if;
  if (x1 ne x2) and (y1 eq y2) then 
    h:=yP - y1;
  end if;
  if IsZero(Q1) eq true then 
    h:=xP - x2;
  end if;  
  if IsZero(Q2) eq true then 
    h:=xP - x1;
  end if;  
  if (IsZero(Q1) eq false) and (IsZero(Q2) eq false) and (x1 ne x2) and (y1 ne y2) then 
    l:=(y1 - y2)/(x1 - x2); 
    s:=l^2;
    h:=(yP - y1 - l*(xP - x1))/(xP + x1 + x2 - s);
  end if;  
  if (Q1 eq Q2) then 
    l:=(3*x1^2 + a) / (2*y1); s:=l^2;
    h:=(yP - y1 - l * (xP - x1))/(xP + x1 + x1 - s);
  end if;  
  return h;
end function;

//Miller's function point doubling: h_{Q,Q}(P) 
hDBL:=function(Q, P, Eda, Edd)
  X1:=Q[1]; Y1:=Q[2]; Z1:=Q[3]; T1:=Q[4];
  A:=X1^2;
  B:=Y1^2;
  C:=Z1^2;
  D:=(X1 + Y1)^2;
  E:=(Y1 + Z1)^2;
  F:=D - (A + B);
  G:=E - (B + C);
  H:=Eda*A;
  I:=H + B;
  J:=C - I;
  K:=J + C;
  X3:=F*K;
  Y3:=I*(B - H);
  Z3:=I*K;
  T3:=F*(B - H);
  CZ2:=2*Y1*(T1 - X1);
  CXY:=2*J + G;
  CXZ:=2*(Eda*X1*T1 - B);
  R:=Vector(4,[X3,Y3,Z3,T3]); 
  h:=CZ2*P[1] + CXY*P[2] + CXZ;
  return R, h;
end function;

//Miller's function point addition h_{S,Q}(P) 
hADD:=function(S, Q, P, Eda, Edd)
  X1:=S[1]; Y1:=S[2]; Z1:=S[3]; T1:=S[4];
  X2:=Q[1]; Y2:=Q[2]; Z2:=Q[3]; T2:=Q[4];
  A:=X1*X2;
  B:=Y1*Y2;
  C:=Z1*T2;
  D:=Z2*T1;
  E:=D + C;
  F:=(X1 - Y1)*(X2 + Y2) + B - A;
  G:=B + Eda*A;
  H:=D - C;
  I:=T1*T2;
  X3:=E*F;
  Y3:=G*H;
  Z3:=F*G;
  T3:=E*H;
  CZ2:=(T1 - X1)*(T2 + X2) - I + A;
  CXY:=X1*Z2 - X2*Z1 + F;
  CXZ:=(Y1 - T1)*(Y2 + T2) - B + I - H;
  R:=Vector(4,[X3,Y3,Z3,T3]); 
  h:=CZ2*P[1] + CXY*P[2] + CXZ;
  return R, h;
end function;

//final exponentiation k = 8, D = 3
finalexp:=function(f, x0, p, q, Fp)
  t0:=1;
  f:=Frobenius(f, Fp, 4)/f; 
  x2:=x0^2;
  x3:=x2^2;
  x4:=x3^2;
  x5:=x0*x2;
  l:=(3*x4 - 6*x3 + x2 - 2*x0 + 7) div 3;
  y:=f^(l);
  t:=Frobenius(y, Fp, 3);
  t0:=t0*t;
  y0:=y^(x0*(x3 - 1))*f^x2 ;
  t:=Frobenius(y0, Fp, 2);
  t0:=t0*t;
  y1:=f*y^x5;
  y2:=Frobenius(y1, Fp, 4);
  t:=Frobenius((y2)^x5, Fp, 1);
  t0:=t0*t*y1;
  return t0;
end function; 

//Optimal ate pairing calculation
OptAtePairing:=function(Q, P, QQ, V, H, q, d, m, p, k, z0, Fp)
  f:=1;  
  T:=Q;  
  s:=Intseq(Abs(V[2]),2); 
  for i:=#s - 2 to 0 by -1 do
    T,h:=hDBL(T, P, d, m); 
    f:=h*f^2; 
    if s[i+1] eq 1 then
      T,h:=hADD(T, Q, P, d, m); 
      f:=h*f; 
    end if;
  end for;
  F:=1;  
  T:=QQ;  
  for i:=#s - 2 to 0 by -1 do
    T,h:=hDBL(T, P, d, m); 
    F:=h*F^2; 
    if s[i+1] eq 1 then
      T,h:=hADD(T, QQ, P, d, m); 
      F:=h*F; 
    end if;
  end for;
  f:=f^(V[2])*Frobenius(f, Fp, 1)*F*H;
  f:=finalexp(f, z0, p, q, Fp);
  return f;
end function;

//MAIN PROGRAM
t2:= 0 ; LOOP:= 80 ;
D := 3 ; k := 8 ;                 
p := 13425986401207647338085857063049211160539773205065792188361162406275649328050239575129752958176629396732162268662817430936913029653257717133386936212850769 ;
q := 115870558819778060730706881586134551090987488334851728304740723586960557890161 ;      
h := 115870558819778060730706881586134551090647090686769697543000867243293884681648 ;    
z0:= 4295331013 ; 
a:= 0 ; 
b:= 1 ;
V:=Vector(4,[z0^2,z0,1,0]);

s:=k div 2;
Fp := GF(p);
P<x> := PolynomialRing(Fp);

E:=EllipticCurve([Fp!(a), Fp!(b)]); 
f:=x^3 + a*x + b; 
n:=h*q; 
t:=p + 1 - n; 
pow:=(p^k - 1) div q;

ROOTS:=Roots(f);
ALPHA:=Vector(3,[ ROOTS[1][1],ROOTS[2][1],ROOTS[3][1] ]);

for i:= 1 to 3 do
  test,sqrt:=IsSquare(1/(3*ALPHA[i]^2 + a));
  if (test eq true) then
    ss:=sqrt;
    alpha:=ALPHA[i];  
    break i;
  end if;
end for;

//Montgomery Curve
A:=3*alpha*ss; 
B:=ss;
//Twisted Edwards Curve
Eda:=Fp!((A + 2)/B); 
Edd:=Fp!((A - 2)/B);

P<z> := PolynomialRing(Fp);
g:=IrreduciblePolynomial(Fp,s); 
Fps<u> := ExtensionField<Fp,z|g>;
R<zz> := PolynomialRing(Fps);

for i1:= 1 to 5 do
  for i0:= 0 to 5 do
    fX:=i1*u + i0; 
    FX:=zz^(k div s) - fX;
    if IsIrreducible(FX) eq true then
      Fpk<v> := ExtensionField<Fps,zz|FX>; 
      x0:=Roots(FX, Fpk); 
      w:=x0[1][1];
      break i1;
    end if;
  end for;
end for;

Et:=EllipticCurve([Fps!(a*w^4),Fps!(b*w^6)]);
Zt<yy> := PolynomialRing(Integers());
Nt := Resultant(yy^k - 1, yy^2 - t*yy + p);
ht := Nt div q^2;

alphat:=Fps!(alpha*w^2); 
sst:=Fps!(ss/w^2);
//Montgomery Curve
At:=3*alphat*sst; 
Bt:=sst;
//Twisted Edwards Curve
Edat:=Fps!((At + 2)/Bt); 
Eddt:=Fps!((At - 2)/Bt);

Ek:=EllipticCurve([Fpk!(a), Fpk!(b)]); 
Zk<xx> := PolynomialRing(Integers());
Nk := Resultant(xx^k - 1, xx^2 - t*xx + p);
hk := Nk div q^2;

Edak:=Fpk!(Eda); 
Eddk:=Fpk!(Edd); 
alphak:=Fpk!(alpha); 
ssk:=Fpk!(ss); 

for i:= 1 to LOOP do
//create two points PW and Pa in E(Fp) of order q (two points because of the bilinearity check)
G:=Random(E); 
nP:=h; 
P:=nP*G; 
nP:=Random(1, q - 1); 
Pa:=P; 
PW:=nP*P; 

MP:=MontPoint(PW, alpha, ss); 
EP:=EdPoint(MP);

MPa:=MontPoint(Pa, alpha, ss); 
EPa:=EdPoint(MPa);

eta:=w*(EP[3] + EP[2])/(EP[1]);
y0:=(EP[2]/EP[3]);
PT:=Vector(2,[eta,y0]);

etaa:=w*(EPa[3] + EPa[2])/(EPa[1]);
y0a:=(EPa[2]/EPa[3]);
PTa:=Vector(2,[etaa,y0a]);

//create two points QW and Qa in Et(Fp^4) of order q (two points because of the bilinearity check)
G:=Random(Et); 
nQ:= ht ; 
Q:=nQ*G;  
QW:=Q; 
xQ:=QW[1]*w^(-2); 
yQ:=QW[2]*w^(-3); 
E1:=Ek![xQ,yQ]; 
XQ:=Frobenius(xQ, Fp, 1); 
YQ:=Frobenius(yQ, Fp, 1); 
D1:=Ek![XQ,YQ]; 
QW:=D1-E1; 

xQ:=QW[1]*w^(2);
yQ:=QW[2]*w^(3);
QW:=Et![xQ,yQ]; 

nQ:=Random(1, q - 1); 
Qa:=QW; 
QW:=nQ*QW; 

MQ:=MontPoint(QW, alphat, sst); 
EQ:=EdPoint(MQ); 

MQa:=MontPoint(Qa, alphat, sst); 
EQa:=EdPoint(MQa); 

//QQ = [z0]QW and QQa=[z0]Qa in Et(Fp^4) of order q (for computing f_{x,[x]Q}(P))
QQ:=z0*QW;
QQa:=z0*Qa;

MQQ:=MontPoint(QQ, alphat, sst); 
EQQ:=EdPoint(MQQ); 

MQQa:=MontPoint(QQa, alphat, sst); 
EQQa:=EdPoint(MQQa); 

//Qh and Qha in Ek(Fp^k) of order q (for computing h_{[p^2]Q,[xp]Q}(P))
Qh:=Ek![QW[1]/w^2,QW[2]/w^3]; 
Qha:=Ek![Qa[1]/w^2,Qa[2]/w^3];

Q1:=(V[3]*p^2 + V[4]*p^3)*Qh; 
Q2:=(V[2]*p)*Qh; 
h2:=hh(Q1, Q2, PW, a); 

Q1a:=(V[3]*p^2 + V[4]*p^3)*Qha; 
Q2a:=(V[2]*p)*Qha; 
h2a:=hh(Q1a, Q2a, Pa, a);

H:=h2; 
Ha:=h2a; 

nQP:=nQ*nP;
t1 := Cputime();
opt:=OptAtePairing(EQ, PT, EQQ, V, H, q, Edat, Eddt, p, k, z0, Fp); 
t2:=t2 + Cputime(t1); 
end for;
t2:=t2/LOOP;
printf"average time for %o pairings = %o\n", LOOP, t2;

//Bilinearity check
OptAtePairing(EQ, PT, EQQ, V, H, q, Edat, Eddt, p, k, z0, Fp) eq OptAtePairing(EQa, PTa, EQQa, V, Ha, q, Edat, Eddt, p, k, z0, Fp)^(nQP);